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Abstract. The two-point angular correlation function is a traditional method used 
to search for deviations from expectations of isotropy. In this paper we develop and 
explore a statistically descriptive three-point method with the intended application 
being the search for deviations from isotropy in the highest energy cosmic rays. 
We compare the sensitivity of a two-point method and a "shape-strength" method 
for a variety of Monte-Carlo simulated anisotropic signals. Studies are done with 
anisotropic source signals diluted by an isotropic background. Type I and II errors 
for rejecting the hypothesis of isotropic cosmic ray arrival directions are evaluated for 
four different event sample sizes: 27, 40, 60 and 80 events, consistent with near term 
data expectations from the Pierre Auger Observatory. In all cases the ability to reject 
the isotropic hypothesis improves with event size and with the fraction of anisotropic 
signal. While ~ 40 event data sets should be sufficient for reliable identification of 
anisotropy in cases of rather extreme (highly anisotropic) data, much larger data sets 
are suggested for reliable identification of more subtle anisotropics. The shape-strength 
method consistently performs better than the two point method and can be easily 
adapted to an arbitrary experimental exposure on the celestial sphere. 
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1. Introduction 

Cosmic rays with energies above 10 EeV (10^^ eV) have been observed[l, 2, 3, 4]. 
However, the sources of these cosmic rays (CR) are unknown and the physics responsible 
for accelerating CR to these energies is at best conjecture. Evidence supporting an extra- 
galactic origin of these CR is the observation of energy flux suppression consistent with 
the GZK-effect[5, 6] by the High Resolution Fly's Eye Experiment [3] and the Pierre 
Auger Observatory (Auger) [4]. The primary evidence supporting the astrophysical 
origin of these CR (as opposed to, say, heavy relic decay) is the lack of an observable 
flux of photons by Auger[7, 8, 9] and the lack of neutrinos observed by ANITA[10, 11]. 

If the sources are astrophysical, expectations for asymmetries in the arrival 
directions increase at the very highest CR energies because the local (< 100 Mpc) 
universe is very anisotropic[12, 13] and the GZK-effect[3, 4] at these energies implies 
that the sources are local. Observation of an anisotropy in the arrival directions of CR 
would be an important step towards identifying the sources of these ultra high energy 
particles. 

Evidence for structure (anisotropy) in the arrival directions has been reported 
[14, 15, 16, 17, 18, 19, 20, 21]. The most compelling observational evidence consistent 
with astrophysical expectations of anisotropy is arguably the 27 events with energy 
greater than 57 EeV recently reported by Auger in [18, 19]. Using the Veron-Cetty 
- Veron (VCV) catalog[12], the active galactic nuclei (AGN) maximum redshift and 
correlation angle chosen by Auger deflned a limited area (eflFectively 21%) of the sky [19]. 
Reported at the 1% signiflcance level, the Auger AGN to CR correlation signal is 
evidence for a flux of CR enhanced near known low-redshift extra-galactic objects[19]. 

As even the largest experiments accumulate the very highest energy CR only slowly, 
:j: the development of new, more sensitive, techniques to search for deviations from 
isotropy is of particular interest. In contrast to the catalog dependent method used by 
[18, 19, 22], in this paper we study the effectiveness of two catalog independent methods. 
Catalog independent techniques avoid the penalty factors for scans over many different 
catalogs and/or the need to restrict the CR data based on limited sky coverage of a 
catalog. The flrst catalog independent technique is a binned two-point {2-Pt) angular 
correlation method (section 2.2). We also introduce a new three-point method which 
uses a shape and a strength parameter [S-S^ or Shape-Strength) for each triplet of events 
(section 2.3). Both methods are compared throughout via the binned-likelihood analysis 
described in section 2.1. 

Arguably, the primary impediment to deflnitive CR source identiflcation is the 
small number of ultra-high energy events (those near and above the GZK cut-off, which 
are most likely to be anisotropic). While lower energy events are more abundant, their 
sources are likely to be further away, where the universe is isotropic. Furthermore, 
galactic/intergalactic magnetic flelds are likely to wash out any correlation with 

:|: For example, the Auger event rate for CR above the GZK knee, ~ 56 EeV, is on the order of two 
per mo nth [4]. 
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the sources of lower energy events (neglecting the possible effects of magnetic field 
caustics[23, 24]). Thus, as one decreases the minimum observed energy one expects to 
include events which dilute any high energy anisotropy signal. Furthermore, there is 
typically significant error in the value of an observed energy (as much as 25% [4]). We 
therefore pay careful attention to the performance of the methods under variation in the 
total number of events and dilution factor (signal to isotropic background) for different 
types of signals in section 3.1. 

2. Methods 

The 2-Pt (section 2.2) and S-S (section 2.3) methods are compared using the analysis 
paradigm described in section 2.1 When needed for a concrete example, we use the 
largest currently operating observatory (Auger) for representative data set sizes and sky 
exposure[18, 19]. The methods presented here, however, can be applied to a spherical 
data set of any size and with an arbitrary experimental exposure. 

2.1. Analysis Paradigm 

We use a similar analysis paradigm for both the 2-Pt and S-S methods to calculate 
a p-value for rejecting the isotropic (null) hypothesis, H^g^- Each method uses binned 
parameters to compute a pseudo-log-likelihood test statistic Sp, "pseudo-" because the 
bins are correlated. The correlation does not effect the final answer because the value 
is derived by comparing the distribution of the Hp in a test sky to that of identically 
analyzed isotropic skies. The flatness of the distribution of p- values for isotropic test 
skies has been verified. The parameter space for the 2-Pt method is the angular distance 
between two events. For the S-S method the parameter space is two dimensional. In 
neither case is the parameter space scanned to determine an optimal value. Instead, we 
compare the entire observed distribution to that expected by isotropy. 

For a given set of cosmic ray events (referred to here as a sky) we compute 
Sp by comparing the binned distribution of the test sky's parameter (s) to the 
parameter (s) distribution expected from an isotropic sky. The probability for observing 
Hobs doublets (2-Pt) or triplets (S-S) from the test sky in the i^^ parameter bin, 
given that you expect [25] to see n^xp from an isotropic sky, is approximated [26] by 
a Poisson distribution Pi{nobs\'i^exp) = ^exp^^^^e^^^^^/ngxp!- The pseudo-log-likelihood is 
Sp = ^^lnP^(no6s|nea;p), where the sum is carried out over the bins of the parameter 
space. The ratio of the number of isotropic skies with Sp less than that of the test sky 
to the total number of simulated isotropic skies gives the value for the test sky. 

In the following discussion fk is defined as the arrival direction of the k^^ event in 
a sky. This (unit) vector has Cartesian coordinates {r^;, r^, r^} when projected from the 
galactic sphere. 
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2.2. Two-Point Correlation 

The 2-Pt correlation distribution is calculated by computing the number of event pairs in 
a test sky as a function of the angular distance between any two events, 9 = co8~^ (rj-rk) 
(see [20, 21, 27] for similar methods). We use 5° bins for 9 G [0°,180°), so that the 
pseudo-log-likelihood is the sum over all angular scales, Sp = In Po^Uobsl'^exp) • 

2.3. Shape- Strength 

This method involves an eigenvector decomposition, or principle component analysis, of 
the arrival directions of all sets of triplets found in the data set. It is inspired primarily 
by Fisher [28] (see also [29, 30]) but differs in that we decompose all subsets of triplets 
in a sky to obtain a test statistic. 

For each triplet we calculate the components of the symmetric (3 x 3) orientation 
matrix T[28]. In Cartesian coordinates, Tij = | "^j^^t^ipi^t {^i^j)k ^ 
largest eigenvalue of T, ri, results from a rotation of the triplet about the principle 
axis Ui. The middle and smallest eigenvalues correspond to the major U2 and minor 
Us axis respectively. The left panel of Figure 1 shows a graphical illustration of these 
eigenvectors. The eigenvalues satisfy ti + T2 + = I and ri > r2 > ra > 0, and thus 
there are only two independently varying parameters for any triplet. 

It is convenient and statistically descriptive to work with a shape^ 7, and a strength^ 
parameter [28]; 



As C increases from to oc the events in the triplet become more concentrated. 
Generally, as 7 increases from — oc to +oc the shape of the triplet transforms from 
elliptical, i.e. strings, to symmetric about Ui^ i.e. point source. See the right panel of 
Figure 1 for a schematic representation. In Figure 2 we show the how the variation of 
the ellipticity of a source on the galactic sphere effects the shape-strength parameter 
space. 

To compute the test statistic Sp using this method we sum over sixty bins for 
7 G [—3.0,3.0) and seventy-five bins for ^ G [0.0, 15.0), i.e. Sp = ^^(^\nP^c,{nohs\'^exp)- 
We have checked that this parameter range is sufficient to cover event sets like those 
expected by Auger and that little is gained by enlarging the range. 

3. Results 

In order to gain confidence in, and intuition about, the S-S method we apply it (and 
the 2-Pt correlation) to an astro-physically motivated simulated (mock) data set in 
section 3.1. The results of applying the S-S and 2-Pt methods to the 27 most energetic 
events from Auger[18, 19] are reported in section 3.2. 
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Figure 1. Left: The eigenvectors of a triplet of events on the sphere (5'^) are the 
principle axis ui^ the major axis U2 (pointing into the page) and the minor axis us. 
The eigenvalues of these vectors are used to compute this triplet's shape and strength. 
Right: An intuitive interpretation of the shape and strength parameters. As the 
strength parameter ( increases from to oo, the events become more concentrated. As 
the shape parameter 7 increases form —00 to +00, the events become more rotationally 
symmetric or less elongated. 



3.1. Mock Signals 

In weighing the effectiveness of a method for rejecting the isotropy hypothesis Hiso for 
a given CR sky we are interested in the probabihties for two types of testing errors[31]. 
A type I error is the probabihty a (significance or p-value) of rejecting H^so given that 
Hiso is true; in practice it should be chosen a priori. In this analysis we choose the 1% 
significance level. One percent is arbitrary and is chosen to be the same as the value 
used in [18, 19]. One could choose, for example, 0.1% but this would require more 
data and/or a higher fraction of anisotropic signal to be detected. For each method 
this choice corresponds to a unique Sp^; we determine Sp^ such that the ratio of the 
number of isotropic skies with Sp less than Sp^ is a = 1%. We use the 10^ isotropic 
skies to determine the upper bound of the signal region of likelihood space. 

A type II error is the probability, /?, of accepting H^so (i.e. of rejecting the mock, or 
toy, signal hypothesis ^sig) given that H^s^ is false. This value is dependent on the choice 
of ^sig- As we are interested in the effectiveness of accepting the signal hypothesis, we 
use the quantity 1 — /?, called the power of the method[31]. By applying the ensemble 
of each mock signal to each method we estimate the power as the ratio of mock signal 
skies with Sp < Sp^ to the total number of mock skies. As a heuristic measure we 
will describe a method's power as "good" if it is at least 90%, i.e. a high probability 
(1 — /? > 0.90) to observe an anisotropy when there is an anisotropy in the data, and 
questionable if it is less than 90%. 

Of significant physical interest is the ability of these catalog independent methods 
to detect signals generated from a catalog dependent map. To this end, we have studied 
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simulated data sets generated from subsets of the VCV[12] galaxy catalog. We consider 
only galaxies with reshift z < 0.02 and we weight each galaxy either by a acceptance 
factor or not at all. We simulate events arriving from these galaxies with a random 
component given by a two dimensional Gaussian centered on the galaxy and with 
deviation a = 3°. These choices are arbitrary in the sense that they describe some 
subset of nearby AGN with events smeared by a few times the angular resolution of 
Auger [19]. It should also be noted that the redshift weighted map (see Figure 3) is 
highly anisotropic, consisting of a number of small to medium scale clumps on the 
celestial sphere, and is likely to yield multiple events per sky within these groups. In 
contrast, the unweighted VCV map (see Figure 4) is notably more dispersed on the 
sphere. 

The true CR data is likely to contain a mixture of background events and signal 
events. To explicitly study this dilution effect has on the power we separately construct 
mock ensembles in which each sky has a certain ratio, r, of signal events to the total 
number of events, with r = 0.2, 0.4, . . . , 1.0. Notice that, because our methods use all 
the triplets or doublets in a given sky, the mixture Sp distributions are not a simple 
sum of the signal and isotropic Hp distributions. 

Detection power is also strongly effected by the number of (high energy) events in a 
sky. The effect can be similar to those of signal dilution in that the power is decreased. 
We generate ensembles of 10^ skies with 27, 40, 60, and 80 events per sky from the 
VCV catalog. Results for all combinations of source purity and number of (mock) data 
events are reported in Figures 3 and 4. The dark blue regions in lower plots of Figure 3 
show that at least 40 — 80 events with (60 — 40)% signal is required to achieve a high 
detection power, 1 — /? ~ 90%, for the redshift weighted VCV maps. The un- weighted 
VCV maps in Figure 4 require a nearly pure signal and 60 or more events to have high 
detection power. 

In general, where one method is good (power, 1 — /? > 90%) so is another; the 
methods are correlated. However, the S-S method consistently performs better than the 
two point correlation for the types of signals discussed here. 

3.2. Auger Data 

It is of interest to apply these techniques to experimental data. The largest public ultra 
high energy data set is the 27 events that form the basis of the Auger result reporting 
evidence for anisotropy (at the 1% significance level) in the highest energy CRs[18, 19]. 
The p-values obtained are: p ^ 3% for the 2-point method§ and p ^ 0.2% for the S-S 
method. Thus, of these two methods only the S-S method would pass a requirement of 
p < 1% as evidence of anisotropy. Note that these events are known to be anisotropic 
- by the methods described in [18, 19, 21] - and therefore the p-values reported here 
reflect only on the statistical methods described in this paper. 

§ We note that the 2-point method used here differs from the auto-correlation analysis performed on 
the 27 events in [21]. 
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4. Conclusion 

In this paper we have introduce a shape-strength method for testing isotropy on the 
unit sphere. We have shown that this method uses pattern-descriptive parameters and 
can consistently out-perform the two-point correlation method. By simulating artificial 
and astrophysically motivated signals of various sizes and purity we can gauge how 
this method might perform on real data. The S-S method out performs the two-point 
method for all of our parameter choices. 

The S-S method was found to detect the redshift weighted VCV toy signal (having 
significant small scale anisotropics) with at least ~ 50% signal purity and about 60 
events in > 90% of the Monte Carlo simulations. We also wish to emphasize from the 
analysis of the diluted mock signals that when the signal to all ratio r > 50% we can 
expect that a redshift weighted "VCV-like" CR signal should be identified with power 
> 50% by both methods for data sets with > 40 events. The unweighted VCV toy 
signal (which is more diflFuse on the sphere) is only reliably detected with greater than 
80 events and 80% signal purity. 

In agreement with qualitative expectations, this analysis demonstrates quantita- 
tively how both signal purity and the total number of events dramatically effect the 
signal detection power. Furthermore, while sources types with significant small scale 
anisotropy can be detected with modest signal purity and total number of events, anal- 
ysis of more subtle anisotropics suggest that either high purity signals or, more likely, 
much more data are needed for reliable identification. 
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Figure 2. Le/t column: Histogram of 10^ skies of 27 Monte-Carlo cosmic rays 
simulated from a single source centered on {/, b} = {—30.0, 0.0} in galactic coordinates. 
We use the Fisher-Bingham distribution [33] on the sphere with = 400.0 to generate 
these events. For the spherically symmetric point-like (top) distribution we use = 0.0. 
For the elliptically shaped (bottom) distribution we use /3 = 200.0 with the major axis 
pointing perpendicular to the galactic plane. See [33] for a detailed description of the 
parameters k and f3. Right column: The ensemble average (over all 10^ sets of 27 event 
skies) of the \nP{nobs\^exp) parameter space of the S-S method for the point-like (top) 
and elliptically shaped (bottom) toy anisotropics. In the bottom (right) panel one 
can see the relatively small deficit of triplets generated from the source with 7^1 
and ( ^ 2 in addition to the large excess of triplets with 7^0 and ( ^ 8. The 
deficit arises from the non-uniform isotropic exposure of Auger [18, 19] and the excess 
from the simulated source. Both features contribute to the pseudo-likelihood where 
no distinction is made between excess and deficit of triplets. These two signals can be 
consistently detected with both the 2-Pt and S-S methods. 
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Figure 3. Top left: Histogram of 10^ skies of 27 Monte-Carlo cosmic rays simulated 
from the VCV[12] catalog. We select objects with redshift Zmax < 0.020 and they 
are weighted by Each simulated CR is drawn from a collection of 2D-Gaussian 

probability distributions centered on the the catalog sources, with deviation a = 3°. 
(See section 3.1.) Top right: Using the VCV catalog we plot the ensemble average of 
the liiP{nobs\'i^exp) parameter space of the S-S method. Bottom row: Using the VCV 
ensemble files we can study the detection power 1 — ^ as a function of the number of 
events per sky and fraction of each sky containing signal events using both the 2-Pt 
(left) and the S-S {right) methods. 
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Figure 4. Top left: Histogram of 10^ skies of 27 Monte- Carlo cosmic rays simulated 
from the VCV[12] catalog. We select objects with redshift Zmax < 0.020 and they 
are not weighted. Each simulated CR is drawn from a collection of 2D-Gaussian 
probability distributions centered on the the catalog sources, with deviation a = 3°. 
(See section 3.1.) Top right: Using the VCV catalog we plot the ensemble average of 
the liiP{nobs\'i^exp) parameter space of the S-S method. Bottom row: Using the VCV 
ensemble files we can study the detection power 1 — ^ as a function of the number of 
events per sky and fraction of each sky containing signal events using both the 2-Pt 
(left) and the S-S {right) methods. 
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